Understanding 911 calls around UW

This analysis adapte from Walking the Beat: Mining Seattle’s Police Report Data by Jeff Wong

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

1. Including a few packages

if (!require(data.table)) { install.packages('data.table'); require(data.table) }

Loading required package: data.table
data.table 1.9.6  For help type ?data.table or https://github.com/Rdatatable/data.table/wiki
The fastest way to learn (by data.table authors): https://www.datacamp.com/courses/data-analysis-the-data-table-way

Attaching package: ‘data.table’

The following object is masked _by_ ‘.GlobalEnv’:

    .N

if (!require(ggplot2)) { install.packages('ggplot2'); require(ggplot2) }

Loading required package: ggplot2
package ‘ggplot2’ was built under R version 3.3.2

if (!require(reshape2)) { install.packages('reshape2'); require(reshape2) }

Loading required package: reshape2

Attaching package: ‘reshape2’

The following objects are masked from ‘package:data.table’:

    dcast, melt

if (!require(devtools)) { install.packages('devtools'); require(devtools) }

Loading required package: devtools

if (!require(ggmap)) { install.packages('ggmap'); require(ggmap) }

Loading required package: ggmap
Google Maps API Terms of Service: http://developers.google.com/maps/terms.
Please cite ggmap if you use it: see citation('ggmap') for details.

Attaching package: ‘ggmap’

The following object is masked _by_ ‘.GlobalEnv’:

    crime

if (!require(glmnet)) { install.packages('glmnet'); require(glmnet) }

Loading required package: glmnet
package ‘glmnet’ was built under R version 3.3.2Loading required package: Matrix
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loaded glmnet 2.0-13

if (!require(fields)) { install.packages('fields'); require(fields) }

Loading required package: fields
package ‘fields’ was built under R version 3.3.2Loading required package: spam
package ‘spam’ was built under R version 3.3.2Loading required package: dotCall64
Loading required package: grid
failed to assign RegisteredNativeSymbol for transpose to transpose since transpose is already defined in the ‘spam’ namespacefailed to assign RegisteredNativeSymbol for toeplitz to toeplitz since toeplitz is already defined in the ‘spam’ namespaceSpam version 2.1-1 (2017-07-02) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: ‘spam’

The following objects are masked from ‘package:base’:

    backsolve, forwardsolve

Loading required package: maps
package ‘maps’ was built under R version 3.3.2

if (!require(dplyr)) { install.packages('dplyr'); require(dplyr) }

Loading required package: dplyr
package ‘dplyr’ was built under R version 3.3.2
Attaching package: ‘dplyr’

The following objects are masked from ‘package:data.table’:

    between, last

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

if (!require(lubridate)) { install.packages('lubridate'); require(lubridate) }

Loading required package: lubridate

Attaching package: ‘lubridate’

The following objects are masked from ‘package:data.table’:

    hour, mday, month, quarter, wday, week, yday, year

The following object is masked from ‘package:base’:

    date

2. Loading the data

MGH_LONGITUDE = -122.307987
MGH_LATITUDE = 47.655038
PATH_DATA_ALL = paste(getwd(), "/data/Seattle_Police_Department_911_Incident_Response.csv", sep = "")
PATH_DATA_MGH = paste(getwd(), "/data/Seattle_Police_Department_911_Incident_Response_Near_MGH.csv", sep = "")
PATH_DATA_MGH_TINY = paste(getwd(), "/data/Seattle_Police_Department_911_Incident_Response_Near_MGH_tiny.csv", sep = "")
if(file.exists(PATH_DATA_MGH)){ # switch variable to PATH_DATA_MGH_TINY if there are issues with scale of data
  data <- fread(PATH_DATA_MGH, header=T, sep=",")
} else {
  data_all = fread(PATH_DATA_ALL, header=T, sep=",")
  data_all <- data_all %>% mutate(dist_mgh = sqrt((MGH_LATITUDE-Latitude)^2 + (MGH_LONGITUDE-Longitude)^2)) # Finding crimes closest to MGH
  data_mgh <- filter(data_all, dist_mgh<0.01) #22,110 rows
  data_mgh_tiny <- filter(data_mgh, dist_mgh<0.007) #5,100 rows
}

C function strtod() returned ERANGE for one or more fields. The first was string input '4.94077829302447e-314'. It was read using (double)strtold() as numeric value 4.9407782930244673E-314 (displayed here using %.16E); loss of accuracy likely occurred. This message is designed to tell you exactly what has been done by fread's C code, so you can search yourself online for many references about double precision accuracy and these specific C functions. You may wish to use colClasses to read the column as character instead and then coerce that column using the Rmpfr package for greater accuracy.

#adding columns to the data
data[,at_scene_time_ts := as.POSIXct(strptime(`At Scene Time`, "%m/%d/%Y %I:%M:%S %p"))] #converting time from String to date and time representation (POSIXct)

data[,at_scene_time_hr := hour(ymd_hms(as.character(at_scene_time_ts)))]

data[,at_scene_time_date := as.Date(at_scene_time_ts)]

data[,at_scene_time_week := floor(as.numeric(at_scene_time_date - min(at_scene_time_date, na.rm=T)) / 7) + 1]

data[,event_clearance_ts := as.POSIXct(strptime(`Event Clearance Date`, "%m/%d/%Y %I:%M:%S %p"))]

data[,event_clearance_date := as.Date(event_clearance_ts)]

data[,event_clearance_hr := hour(ymd_hms(as.character(event_clearance_ts)))]

data[,time_until_event_clear := as.numeric(event_clearance_ts - at_scene_time_ts)]

data[,`Initial Type Group` := factor(`Initial Type Group`)]

data[,`Event Clearance Group` := factor(`Event Clearance Group`)]

data[,`Zone/Beat` := factor(`Zone/Beat`)]

data[,LatitudeBin := round(Latitude, 3)]

data[,LongitudeBin := round(Longitude, 3)]

data[,crime_type := ifelse(`Event Clearance Group` %in% crimes.violent, "Violent", ifelse(`Event Clearance Group` %in% crimes.serious, "Serious", "Minor"))]

View(data)

3. Missing Data

Missing data can be a problem. If a large proportion of data is missing, we may end with results that are not representative of the population.

TODO There are several options for date/time (at_scene_time_ts, event_clearance_ts). Figure out what proportion of reports have values. hint: is.na can be helpful

total_reports = nrow(data)
total_reports

[1] 22110

# write code below!

4. Frequency of Crimes by Day of Week

Let’s do a quick sanity check of our data by looking at the number of reports by day of week. Because police services are available every day of the week, we would expect at least some reports occuring each day of the week. Run the code below to verify that.

sanityCheck = data.frame(data[,table(weekdays(event_clearance_date))]); # frequency of crimes by day of week
View(sanityCheck)
sanityCheck$Var1 = factor(sanityCheck$Var1, levels = rev(c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday")))
ggplot(sanityCheck,
       aes(x = Var1, y = Freq, fill = as.numeric(Var1) %% 2)) +
  geom_bar(stat = 'identity') +
  coord_flip() +
  theme_bw() + 
  xlab ("Day of Week") + ylab("Count") +
  guides(fill=F)

5. Frequency of Crimes by Time of Day

Perhaps we should have more officers patroling during times of day which are high crime. TODO Plot the number of reported incidents by hour of day and see what suggestions you may make.

hint this code should look similar to the previous block (Frequency of Crimes by Day of Week). challenge: Crime by time of day may vary by weekend days and weekdays. Plot the crime by times of day for weekdays and weekend days separately.

# Your code below!

6. Types of contacts

Not all crime is the same! Let’s see the frequency of each type of crime. We provide a function ggplot.freqtable.1d to generate the plot. Run the code below.

#Uses a frequency table to build a bar chart that is ordered by the counts
ggplot.freqtable.1d = function(x) {
  freq = data.frame(x); colnames(freq)[1] = "label"
  freq.filtered = droplevels(subset(freq, label != "" & !is.na(label)))
  freq.filtered$label = factor(freq.filtered$label,
                               levels = as.character(freq.filtered[order(freq.filtered$Freq, decreasing = F),]$label),
                               ordered = F)
  ggplot(freq.filtered,
         aes(x = label, y = Freq, fill = label)) +
    geom_bar(stat = 'identity') +
    coord_flip() +
    theme_grey()
}
#Ploting frequency of each type of report
ggplot.freqtable.1d(data[,table(`Event Clearance Group`)]) +
  xlab("Contact Type") + ylab("Count") +
  theme_bw() +
  guides(fill=F)

7. Digging deeper in the type of event

There are a lot of “Disturbances” and “Suspicious Circumstances” reported, but what does that mean? The “Event Clearance Description” field may provide more information.

TODO: Using the ggplot.freqtable.1d function from above, generate plots which provide more information on the types of crimes which are labeled as “Disturbances.” hint: Filter the data by Event Clearance Group challenge: Do the same analysis for the “Suspicious Circumstances” group.

8. Heatmap of Crimes

Some places may be more of hotbeds from crime. Run the code below to figure out which areas have more reported crimes.

ZOOM <- 15 # number between 1 and 21
seattle = get_map(location = c(lon = MGH_LONGITUDE, lat = MGH_LATITUDE), zoom = ZOOM, maptype = 'roadmap')

Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.655038,-122.307987&zoom=15&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false

data.filtered = data[!is.na(time_until_event_clear) & time_until_event_clear > 0,
                     list(count = .N),
                     list(LongitudeBin, LatitudeBin)]
# If you run into an error about GeomRasterAnn being built on an incompatible version of ggproto, you must reinstall ggmap. Use this command: install.packages("ggmap", type = "source")
ggmap(seattle) +
  geom_point(data = data.filtered[count > 10],
             mapping = aes(x = LongitudeBin, y = LatitudeBin, size = log10(count), alpha = log10(count)), color = 'blue') +
  theme_grey() +
  theme(legend.position="top") +
  xlab("Longitude") + ylab("Latitude") + labs(title = "Heatmap of Crimes")

9: Responding to casualties

Crimes where people are injured require a fast response. Run this code to look at a heatmap of the time it takes to clear crimes involving injury. TODO: Change the type of crime to determine if clearance time varies by crime.

crime = 'PERSON DOWN/INJURY'
data.filtered = data[!is.na(time_until_event_clear) & time_until_event_clear > 0]
data.filtered.contact = data.filtered[`Event Clearance Group` == crime]
time_until_event_clear.q01 = quantile(data.filtered.contact$time_until_event_clear, .01); time_until_event_clear.q99 = quantile(data.filtered.contact$time_until_event_clear, .99)
data.filtered.contact = data.filtered.contact[time_until_event_clear > time_until_event_clear.q01 & time_until_event_clear < time_until_event_clear.q99]
data.filtered.contact2 = data.filtered.contact[,list(Hours = mean(time_until_event_clear / 3600), count = .N), list(Longitude = LongitudeBin, Latitude = LatitudeBin)]
ggmap(seattle) +
  geom_point(data = data.filtered.contact2,
             mapping = aes(x = Longitude, y = Latitude, color = Hours, alpha = Hours, size = count)) +
  scale_colour_gradientn(colours=c("blue", "red")) +
  scale_alpha(range = c(0.2, 0.7)) +
  scale_size_continuous(range = c(5,20)) +
  guides(alpha=F,size=F) +
  theme_grey() +
  theme(legend.position="top") +
  xlab("Lon") + ylab("Lat") + 
  labs(title = sprintf("Time Until %s Cleared", crime))